Take-home_Ex10

Author

Yuheng Liang

Published

October 28, 2024

Modified

October 28, 2024

Take home exe 03

Data

Importing data and package

package

  • spNetwork, which provides functions to perform Spatial Point Patterns Analysis such as kernel density estimation (KDE) and K-function on network. It also can be used to build spatial matrices (‘listw’ objects like in ‘spdep’ package) to conduct any kind of traditional spatial analysis with spatial weights based on reticular distances.
  • sf package provides functions to manage, processing, and manipulate Simple Features, a formal geospatial data standard that specifies a storage and access model of spatial geometries such as points, lines, and polygons.
  • tmap which provides functions for plotting cartographic quality static point patterns maps or interactive maps by using leaflet API.
  • spatstat, which has a wide range of useful functions for point pattern analysis. In this hands-on exercise, it will be used to perform 1st- and 2nd-order spatial point patterns analysis and derive kernel density estimation (KDE) layer.
  • raster which reads, writes, manipulates, analyses and model of gridded spatial data (i.e. raster). In this hands-on exercise, it will be used to convert image output generate by spatstat into raster format.
  • maptools which provides a set of tools for manipulating geographic data. In this hands-on exercise, we mainly use it to convert Spatial objects into ppp format of spatstat.
pacman::p_load(sf, raster, spatstat, tmap, tidyverse,sparr,spNetwork,dplyr,animation)

data

drug_case <- read_csv("data/2016-01-01-2024-06-30-Philippines.csv")%>%
    st_as_sf(coords = c("longitude", "latitude"), crs = 4326) %>%
    st_transform(crs = 32651) %>%
  mutate(event_date = dmy(event_date)) %>%
  mutate(event_month = year*100 + month(event_date)) %>%
  mutate(event_quarter = year*10 + quarter(event_date)) 
Rows: 6921 Columns: 31
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (23): event_id_cnty, event_date, disorder_type, event_type, sub_event_ty...
dbl  (8): year, time_precision, iso, latitude, longitude, geo_precision, fat...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
write_csv(drug_case,"data/drug_case.csv")

import the boundary data

ph_sf = st_read(dsn = "data/phl_adm_psa_namria_20231106_shp",layer = "phl_admbnda_adm2_psa_namria_20231106")
Reading layer `phl_admbnda_adm2_psa_namria_20231106' from data source 
  `/Users/liangyuhang/Downloads/Maaaaaaaaaark/IS415_g/Take-home_Ex/Take-home_Ex03/data/phl_adm_psa_namria_20231106_shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 88 features and 13 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 114.2779 ymin: 4.587294 xmax: 126.605 ymax: 21.12189
Geodetic CRS:  WGS 84

chack the crs

st_crs(ph_sf)
Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

transform the crs

ph_sf = st_transform(ph_sf,crs = 32651)
st_crs(ph_sf)
Coordinate Reference System:
  User input: EPSG:32651 
  wkt:
PROJCRS["WGS 84 / UTM zone 51N",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 51N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",123,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Between 120°E and 126°E, northern hemisphere between equator and 84°N, onshore and offshore. China. Indonesia. Japan. North Korea. Philippines. Russian Federation. South Korea. Taiwan."],
        BBOX[0,120,84,126]],
    ID["EPSG",32651]]
ph_sf_owin <-as.owin(ph_sf)

show the point

tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(drug_case)+tm_dots()
tmap_mode('plot')
tmap mode set to plotting
hist(drug_case$year)

hist(drug_case$event_quarter)

hist(drug_case$event_month)

  • 1 我应该基于那些点进行分析
  • 2 先进行随机点的分析看菲律宾这些点是不是都是随机的
  • 3 如果不是随机的那么哪里是案件发生比较高的地方
  • 4 热点区域在多长时间内保持稳定？是否存在周期性的热点区域转移
  • 5 是否存在事件上的关联（不同类型的事件上是否会发生转变） 分析：
  • 1 时间趋势上的分析
  • 2 不同时间类型在不同时间上的分析
  • 3 空间热点上的分析

Computing Stkde by year

drug_cases_year = drug_case %>% select(year)
drug_cases_year_ppp = as.ppp(drug_cases_year)
any(duplicated(drug_cases_year_ppp))
[1] TRUE
drug_cases_year_ppp_jit <- rjitter(drug_cases_year_ppp, 
                             retry=TRUE, 
                             nsim=1, 
                             drop=TRUE)
any(duplicated(drug_cases_year_ppp_jit))
[1] FALSE
drug_cases_year_owin = drug_cases_year_ppp_jit[ph_sf_owin]
st_kde_year = spattemp.density(drug_cases_year_owin)
Calculating trivariate smooth...Done.
Edge-correcting...Done.
Conditioning on time...Done.
#plot using 2 for loops, then pump into a variable to hold it, then use animation? Just plot(st_kde) will plot 20165, which doesnt exist
yrs = c(2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024)

animation::saveGIF(
  for(y in yrs){
      plot(st_kde_year,
           main = paste("Drug Cases on ", y)) 
  },
  movie.name = "drug_stkde_year.gif", interval = 0.1, ani.width = 600
)
Output at: drug_stkde_year.gif
[1] TRUE

Computing Stkde by Month

drug_case_month = drug_case %>% 
  select(event_month)
drug_case_month_ppp = as.ppp(drug_case_month)
any(duplicated(drug_case_month_ppp))
[1] TRUE
drug_case_month_ppp_jit <- rjitter(drug_case_month_ppp, 
                             retry=TRUE, 
                             nsim=1, 
                             drop=TRUE)
any(duplicated(drug_case_month_ppp_jit))
[1] FALSE
drug_case_month_owin = drug_case_month_ppp_jit[ph_sf_owin]
st_kde_month = spattemp.density(drug_case_month_owin)
Calculating trivariate smooth...Done.
Edge-correcting...Done.
Conditioning on time...Done.
yrs = c(2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024)
mths = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 ,12)

animation::saveGIF(
  for(y in yrs){
    if(y == 2024){
      mths = c(1, 2, 3, 4, 5, 6)
    }
    for(m in mths){
      plot(st_kde_month, y*100+m,
           main = paste("Drug Cases on ", y, ", ", m,"month"))
    }
  },
  movie.name = "drug_stkde_month.gif", interval = 0.1, ani.width = 600
)
Output at: drug_stkde_month.gif
[1] TRUE